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ABSTRACT 

We used unprecedentedly large 2D and 3D hybrid (kinetic ions - fluid electrons) simulations of non- 
relativistic collisionless strong shocks in order to investigate the effects of self-consistently accelerated 
ions on the overall shock dynamics. The current driven by suprathermal particles streaming ahead of 
the shock excites modes transverse to the background magnetic field. The Lorentz force induced by 
these self- amplified fields tends to excavate tubular, underdense, magnetic-field-depleted cavities that 
are advected with the fluid and perturb the shock surface, triggering downstream turbulent motions. 
These motions further amplify the magnetic field, up to factors of 50-100 in knot-like structures. 
Once downstream, the cavities tend to be filled by hot plasma plumes that compress and stretch 
the magnetic fields in elongated filaments; this effect is particularly evident if the shock propagates 
parallel to the background field. Highly-magnetized knots and filaments may provide explanations for 
the rapid X-ray variability observed in RX J 1713. 7-3946 and for the regular pattern of X-ray bright 
stripes detected in Tycho's supernova remnant. 
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1. INTRODUCTION 



Following the pioneering idea of Fermi ( |1949[ ), in 
the late '70s several authors realized that collisionless 
shoc ks are prominent sites for the acceleration of parti 



clesjKrymskii 1977; A xford et aL|1977| |Bell|1978| |Bland-| 



lord & O striker 1978). More recently, the detection of 
narrow X-ray rims m young supernova remnants (SNRs), 
whose origin has been explained as synchrotron emission 
of relativistic electrons radiating in magnetic fields of a 
few hundreds /iG, has provided evidence that the level 
of magnetization at SNR shocks is much larger than in 



the interstellar me dium (see e.g.. Vink fc Laming||2003 
Bamba et aT|20Q5| |Parizot et aL|2UU6p . 



Tins association between particle acceleration and 
magnetic field amplification has been welcomed by the- 
orists for several reasons. First, the super- Alfvenic 
streaming of accelerated particles is predict ed to ex- 
cite s everal different plasma instabilities (e.g., Bell 1978 



2004 ), which can account for the inferred levels of magne 
tization ( at least as an extrapo lation of the linear theory, 
see, e.g., [Caprioli et al.| [2009). Second, the interstellar 
turbulence cannot scatter particles efficiently enough to 
achieve the highest energies measured in Galactic c osmic 



rays, while self-generated magnetic fields can (e.g. Blasi 



et al. 2007). Finally, amplified magnetic fields may rep- 
resent a key ingredient for explaining t he steep spectr a 
inferred from recent 7-ray observations (Caprioli 2 012[ ). 

The details of such an interplay between accelerated 
particles and magnetic fields have not been completely 
understood yet, mostly because of the difficulty in ac- 
counting for the particle-wave coupling in the fully non- 
linear regime. Nevertheless, collisionless shocks are 
mediated by electromagnetic interactions only, so they 
can be modeled by iteratively moving particles on a 
grid according to the Lorentz force and adjusting the 
electromagnetic configuration via Maxwell equations. 
Such a particle-in-cell (PIC) approach provides great 
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Figure 1. Evolution of the plasma density for a 2D, par- 
allel collisionless shock with sonic and Alfvenic Mach numbers 
M s = Ma = 30, as seen in the downstream frame. The total 
computational box measures 10 4 x 10 3 (c/w p ) 2 , only a portion of 
which is showed here. For comparison, the nominal gyroradius of 
an ion with velocity v s h is tl = v s h/w c = SOc/ujp in the upstream 
magnetic field, while the most energetic particles at the end of the 
simulation have gyroradii as large as r^(^ macc ) ~ 300c/oj p . 

insig ht into the properties of collisionless shocks (see, 
e.g., [Stroman et al. 2009; |Ohira et al.||2QQ9| [Riquelme 
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2009, [20101 |Sironi fc Spitkovsky| [20lT 



fc Spitkovskyl 

Niemiec et al. 2012). However, PIC codes are com- 
putationally quite expensive, in that they require the 
plasma and gyration scales of both ions and electrons 
to be resolved. To partially mitigate this problem, one 
can adopt the so-called hybrid approach, in which the 
ions, which drive the dynamics, are treated kinetically 
while the (massless) electrons are treated as a neutral- 
izing fluid. This approach neglects the small electron 
scales and allows the investigation of more macroscopic 
phenomena in much larger (in physical units) computa- 
tional boxes. In particular, it has been widely exploited 
to study ion acceleration at collisionless s hocks in differ- 
ent astrophysical environments (see, e.g., Win ske 1985 



Log[E f(E)] @t=500(B _1 



Lipatov 2002; Giacalo ne et al.|| 19971 Giacalone ^TEllison 
2000| |Giacalone||2004| |Gargat(§ Sp itkovsk7 poT2j )~ 

Here we show the results ol 2D and 3D hybri d simula- 
tions run wi th the non-relativistic code dHybrid (| Gar gate 
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2007). These simulations allow for unprecedent- 



et al. 

edly large computational boxes, especially in the direc- 
tion transverse to the shock velocity, and can provide 
first-principle insight on the multi-dimensional structure 
of non-relativistic collisionless shocks. We focus our at- 
tention on shocks with high-Mach numbers (M > 30); 
this regime is most relevant for SNR shocks, and so far 
has been scarcely studied because of the large dynamical 
range required to follow both thermal and non-thermal 
particles. 

We investigate, for the first time in a self-consistent 
simulation of a collisionless shock with accelerated parti- 
cles, the formation of peculiar structures in the upstream, 
which are shaped as tubular cavities surrounded by a net 
of dense filaments where the magnetic field is significantly 
amplified. That accelerated particles drive an instabil- 
ity that is filamentary in na t ure has already been pu t 
forwaxdby_ peiT| (120041 [20 05 ] Ij lZirakashvili et al.| fl2008D , 
and |Kogachevskii et al.| ( |2012jj , who have run magneto- 
hydrodynamical simulations with a fixed current imposed 
through a p eriodic box seeded wit h an initial turbulence. 

Recently, |Reville fc Bell| ( |2012[ ) have provided an an- 
alytical derivation ol the cavity growth rate, supporting 
their findings with simulations that couple a fluid treat- 
ment of the background plasma with a kinetic description 
of the energetic ions driving the instability. The insta- 
bility was studied in a 2D slab geometry representing a 
section of the upstream perpendicular to the shock nor- 
mal. 

The simulations presented here differ from those in the 
literature in several respects: 1) the ion current is time- 
dependent and self-consistently generated by shock accel- 
eration, and not estimated by assuming the relativistic 
ions to be isotropic in the shock frame, which translates 
into an anisotropy of order ~ v s /c in the upstream refer- 
ence frame (here v s is the shock speed); 2) we do not need 
to seed the fluid with pre-existing turbulence, since accel- 
erated particles generate it by themselves via streaming 
instability; 3) the more physical setup where filamen- 
tation develops while the fluid is advected towards the 
shock allows us to properly study the evolution and the 
spatial features of cavities and filaments; 4) simulating 
the global shock structure allows us to show how, when 
advected through the shock, cavities and filaments affect 
the nature of the discontinuity and the magnetic field 
topology. 




Figure 2. Top panel: Particle energy spectrum f(E) as a func- 
tion of x in units of E s ^ = mv^ h /2 at time t = SOOo;^ 1 . Bottom 
panel: Time evolution of the particle spectrum for x < 2000c/ oj p , 
as seen in the downstream reference frame. For comparison, 
the Maxwellian distribution corresponding to a temperature T = 
0.85E t h/kB is showed as well (dashed line), with E t h = 3/8E s h 
the temperature given by standard jump conditions. 

In ^2] we outline the main features of our 2D and 3D 
simulations for both parallel and oblique shocks, and dis- 
cuss the mechanisms that lead to the formation of cavi- 
ties and filaments. In 33] we discuss some possible obser- 
vational implications ofsuch filamentation in the context 
of the synchrotron emission detected in young SNRs like 
Tycho and RX J1713.7-3946. 

2. HYBRID SIMULATIONS 

Most of the simulations presented here are for a 2D 
parallel shock (i.e., background field Bo aligned with the 
shock velocity = v s h%) with sonic and Alfvenic Mach 
numbers M s = v s h/c s = Ma — v s h/vA = M = 30 in the 
downstream (simulation) reference frame, corresponding 
to M « 40 in the upstream one. In the 2D simulations, 
we include both in-plane and out-of-plane components of 
the ion momenta and electromagnetic fields. 

In our non-relativistic framework, velocities are nor- 
malized to the Alfven speed va, lengths to the ion skin 
depth c/ujp = va/^c-, and times to cj" 1 , with uj p = 
Y // 47rnoe 2 /m and uj c = eB^/mc being the ion plasma 
and cyclotron frequencies (m is the proton mass). Den- 
sity and magnetic fields are measured in units of the up- 
stream initial values, no and Bq. The ion skin depth is 
resolved with two grid cells and the computational box 
measures (L x x L y ) = 10 4 x 10 3 (c/uj p ) 2 . The time step is 
chosen as At = O.OOIcj" 1 to enhance the energy conser- 
vation with a very small Courant number. The shock is 
generated by introducing a perfectly reflecting wall (left 
boundary in the figures) that initially produces counter- 
streaming particles. This configuration is unstable and 
the system promptly forms a propagating sharp discon- 
tinuity with width of a few gyroradii of the downstream 
thermal particles, behind which the cold upstream flow 
is isotropized. 

The global evolution of the shock obtained with the 
present setup is showed in Figure [I] where the density is 
plotted at different times. The shock propagates to the 
right along x axis. Since M ^> 1, the ratio between the 
plasma density behind and ahead of the shock rapidly 
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Figure 3. Density, parallel (B x ), transverse (B y ,B z ) and total 
(Btot) magnetic field, and Alfven velocity va = Btot / V ^Trmn for 
the 2D simulation in Figurehjat t = bOOou^ 1 . All the quantities are 
normalized to their initial values. A filament with Btot ~ 15 — 20£>o 
and some knots where Btot ~ 20 — 40Bo are clearly visible. 

approaches 4, and the post-shock fluid comes to rest in 
the wall (downstream) frame, dissipating the inflowing 
kinetic energy into heat. 

Figure [2] shows a snapshot of the ion distribution 
function at t = 500CJ" 1 and the time evolution of the 
downstream spectrum. According to first-order Fermi 
mechanism, particles scattered back and forth across the 
shock get accelerated effectively: the ion energy spec- 
trum quickly develops a non-thermal power-law tail that, 
for the chosen parameters, comprises ~ 15% of the total 
ion energy. The details of the ion energy spectrum will 
be discussed in a forthcoming paper. 

The peak of the ion Maxwellian distribution is consis- 
tently shifted to lower temperatures by about the same 
amount with respect to the pure gaseous case (i.e., the 
case with no accelerated particles). At the same time, 
since all the particles are coupled through electromag- 
netic interactions, the pressure in non-thermal ions prop- 
agating upstream slows down the incoming fluid produc- 
ing a so-called precursor (top panel in Figure [2| , a dis- 
tinctive feature of shocks modifie d by the back-reaction 



of accelerated particles (see , e.g., Jones fc Ellisonp991 
Malkov fc O'C. DrurypOQl] for reviews) 
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Figure 4. As in Figure [3] but for an oblique shock with Bo in 
the xy-plane, at an angle v — 20° with respect to v^. 



2.1. Filamentation instability 

The most noticeable feature in Figure [ll is the forma- 
tion of low-density cavities ahead of the shock (n/no ^ 
10 -2 ), whose origin can be explained in the following 
way. Particles accelerated at the shock propagate up- 
stream against the incoming fluid (top panel of Figure 
At first, the magnetic field is the background Bo only; 
therefore, ions do not diffuse in pitch angle, giving rise to 
a coherent current J || B . However, the super- Alfvenic 
streaming of the accelerated particles tends to excite 
transverse magnetic modes via plasma instabilities, even- 
tually seeding the upstream medium with a transverse 
£B. The resulting Lorentz force ex —J x £B pushes 
the plasma (and the frozen-in magnetic field) away from 
the region where the current is stronger, also focusing 
the energetic particles and helping to sustain the insta- 
bility (Bell 2QQ5|). The net result is the formation of low- 
density tunnels rilled with supra-thermal ions, modulated 
by the period of the underlying magnetic perturbation 
(see, e.g., the snapshot at t = SOOo;" 1 in Figure [l]). 

A complete discussion of the most unstable modes ex- 
cited in large Mach number parallel shocks is beyond 
the scope of t his paper. The reader may refer to Amato 
& Blasi (20091) for a kinetic description of resonant and 



non-resonants modes excited in the linear regime and to 
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the simulations by, e .g., Stroman et al. ( 2009 ); Riquelme 
fc Spitkovskyl ( |2009[ ); |Gargate et al.| ( |2010| Tfor a study 
of the fully non-linear regime. Note that Bell's non- 
resonant modes are not particularly relevant here since 
their growth rate is typically shorter than ~ cj" 1 , and on 
the advection time-scale « Lb ox /v s h ^> uo~ x such small- 
wavelength modes have already reached saturation, while 
modes resonant w ith the acceler ated particles are still 
growing (see, e.g., Gargate et al. 2010). A sizable frac- 
tion of the energy in the magnetic turbulence is stored 
on scales as large as lO-lOOc/cjp, comparable with, or 
larger than, the gyroradius ~ v s h/u c = 30c /lj p of ions 
with velocity ~ v s h- 

2.2. Growth of cavities and filaments 



An interestin 
of the cavities. 



question is what de termines the size 



Bell 



2012| ) argued that the 
growth of the cavities is suppressed when their size be- 
comes comparable with the gyroradius of the ions car- 
rying most of the current, since the scattering of these 
particles would suppress the instability. In our setup, the 
growth rate of the filamentation may also be limited by 
advection, while the current may increase with time be- 
cause the acceleration becomes more and more efficient. 

In our simulations cavities grow while being advected 
toward the shock and, by comparing different panels of 
Figure [T] we notice that the typical transverse size of 
the cavities impacting the shock increases with time. 
This is qualitatively consistent with a cavity growth rate 

1/2 

T oc \B±_\^ a c C , where \B±\ is the averaged magnitude of 
transverse component of B, and £ acc is the fraction of 
the to tal pressure in accelerated particles ( Reville & Bell 
[2012] ). In our simulations, £ acc and, in turn, \B±\ in- 
crease with time (see the non-thermal tail in the spectra 
of Figure [2|. At t > 400co>~ 1 the transverse size of the 
biggest cavity measures ~ 300c/cj p , and it is comparable 
with the gyroradius of the ion with the highest energy, 
Emax ~ 100E ah , with E 8h = mv 2 sh /2 (Figureg). 

Lower-energy particles resonate with the self^generated 
magnetic turbulence and are effectively scattered, as 
shown in the top panel of Figure [2] The density of ions 
with E < lOi^h drops with the distance from the shock: 
if these particles were not deflected back, they would es- 
cape the box on a timescale shorter than the simulation 
one. By undergoing multiple interactions with the shock, 
ions are efficiently accelerated, so that E max (t) increases 
with time. On the other hand, freshly accelerated ions 
with energies ~ E max do not have enough waves to res- 
onate with and stream more freely, contributing the most 
to the current because of their anisotropy. 

For t > 500CJ" 1 , E max (t) no longer increases due to 
the finite size of the box in the x direction: however, we 
checked that this limitation is removed for larger boxes 
and consistently longer simulation times. Enlarging the 
box in the y direction, instead, does not accommodate 
larger and larger cavities, attesting to the physical origin 
of their diameter. The transverse size of the box must, 
however, be large enough to resolve the gyration of the 
most energetic ions at any time: the smaller the box in 
the y direction, the earlier filamentation is suppressed. 
Previous simulat ions (e.g., 
|Spitkovsky||2012[ ) used smaf 
effect was less evident. 




Giacalone 2004 Gargate 
ei boxes, and, therefore, the 



Figure 5. Snapshot at t = 175ou c 1 of a 3D hybrid simulation of 
a parallel shock with M s = M A = 6 in a 2000 x 200 x 200(c/cc;p) 3 
box (see text for further details). The color code (right colorbar) 
on the box sides shows the particle density in units of no, while 
the vertical slice illustrates a section of the fluid at x = 520c/ oj p , 
i.e., immediately ahead of the shock. In the slice, the grey-scale 
code corresponds to the ion density, while the colored vectors show 
strength and direction of the magnetic field, in units of Bo. Notice 
the correlation between underdensity and low S-field and how the 
magnetic field is mainly along x in the filaments and coiled inside 
the cavities. 



2.3. Magnetic field structure 

In Figure [3] the parallel (B x ), transverse (B y ,B z ) and 
total (Btot) magnetic fields are shown for the 2D case at 
t = bOOuj^ 1 . In particular, the filamentary structure of 
the upstream is very evident in B x , where organized lat- 
eral modulations significantly affect the mean magnetic 
field topology. While at parallel shocks the ion-driven 
streaming instability is expected to produce a perpendic- 
ular component of the field, a change in B x is a non-linear 
effect which is not predicted in the quasi-linear analysis. 
Also, B x tends to be evacuated from cavities and pushed 
into filaments, while the transverse component winds up 
around the cavity itself (Figure [3]). 

There is also another interesting topological effect that 
arises in the scenario presented here. When a cavity 
is advected through the shock, a bubble of relatively 
rarefied and cold (it contains less kinetic energy to be 
dissipated) plasma is produced in the downstream (see, 
e.g., Figure [I] at t = 300CJ" 1 ). Such a configuration is 
Rayleigh-Taylor unstable, and each bubble is promptly 
filled by avlume of hotter and denser plasma, as visible 
in Figure [ll between t = 400 and 500CJ" 1 . These macro- 
scopic motions systematically stretch the magnetic field 
along the cavity axis, as can be seen in B x and B tot in 
Figure [3] 

A similar phenomenology is recovered for oblique 
shocks, too. Figure [4] shows the output of a simulation for 
a shock with an angle of 20° between Bo and v^, and all 
the other parameters fixed as in the parallel case. Cav- 
ities always develop along the direction of the upstream 
background field, being driven by energetic particles spi- 
raling around Bo- However, in this case the shock is less 
corrugated, and downstream cavities are squeezed more 
rapidly, preventing the formation of very long filaments. 
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2.4. 3D simulations 

To better understand the topology of these structures, 
we also ran 3D simulations with dHybrid. Both the shock 
velocity and the computational domain are scaled down 
by a factor of 5 with respect to the 2D case to compensate 
for the higher computational effort without changing the 
ratio of the box size and the gyroradius of typical parti- 
cles, vl oc v s h- More precisely, we set M s = Ma = 6 in a 
2000x200 x 200(c/cOp) 3 box, allowing for the same time 
and space resolution as in the 2D case. 

Figure [5] shows the rendering of the density and field 
structure near the shock in a 3D simulation. The slice 
shows a transverse section of the fluid immediately ahead 
of the shock, where the contrast between filaments and 
cavities may be as large as 5-10 in both the ion density 
(grayscale) and in total magnetic field (color-coded vec- 
tors). There is a clear correlation between regions of en- 
hanced magnetic field and enhanced density, confirming 
the formation of rarefied, low-B cavities surrounded by 
regions of denser plasma permeated by a stronger field. 

An interesting feature clearly visible in 3D simulations 
is that the field is mainly parallel to the shock normal 
in the dense filaments, while inside the cavities the field 
is coiled and develops a transverse component compara- 
ble with, and even l arger than the one along x (see also 
Reville & Bell||2012|). 



3. OBSERVATIONAL CONSEQUENCES 

Very generally, we can assess that filamentation pro- 
vides substantial magnetic field amplification even in the 
pre-shock medium. The magnetic-field strength in some 
filaments may become as large as ~ 5£?o, while cavities 
may become effectively demagnetized. While low-energy 
accelerated ions are focused in the inner regions of the 
cavities, ions with gyroradii comparable with the cavity 
transverse size can feel the effect of the field amplification 
and scatter more efficiently. In particular, the increase 
of B and the decrease of n in the cavities lead to a local 
Alfven velocity significantly larger than the initial one 
(bottom panels of Figures [3] and [4| . This suggests that 
the phase velocity of the magnetic perturbations may be- 
come a non-negligible fraction of the fluid velocity even 
for strong shocks, eventually affecting the scattering and, 
in turn, the spectra of the accelerated particles. In par- 
ticular, the effect may be crucial to expl ain the steep in - 
ferred from 7-ray observations of SNRs ( |Caprioli||2012[ ). 

Also, the thermal plasma is affected by the tilamen- 
tation process: ahead of the shock the temperature be- 
comes several times larger than at the beginning of the 
simulation. Quite interestingly, the pressure in the ther- 
mal plasma and in the magnetic turbulence are almost 
in equipartition along the precursor. 

Due to the non-linear processes developing in the up- 
stream, the magnetic field orientation and the local 
Alfvenic and sonic Mach numbers may vary significantly 
immediately ahead of the shock (Figures [3] and|4| ; differ- 
ent patches of the upstream plasma are therefore shocked 
in different ways, eventually leading to the onset of co- 
herent motions that stretch the field preferentially along 
the cavities, but also turbulent motions that effectively 
stir the downstream fluid on smaller scales. 

The net result is that the post-shock magnetic field 
may become even larger than what is expected from a 



simple compression by a factor ~ 4 of the transverse 
component of the pre-shock field. In Figure[3j in addition 
to an elongated filament where B > 15 — 20£?o, it is 
possible to spot knots with B ~ 30 — 40£?o- Runs with 
larger Mach numbers (M = 50) show that the total B- 
field in knot-like structures can easily reach up to 90- 
lOO^o. These structures may resemble the ones detected 
in RX J1713. 7-3946, where the fields i nferred by the X- 
ray v ariability are as strong as lmG (Uchiyama et al. 
2007), namely a few hundred times stronger than the 



typical interstellar magnetic field. 

In application to SNRs, if the cavity sizes continue to 
grow to the gyroradius of ions with E max , filamentation 
instability might also account for the detec tion of a char- 
acter istic pattern of X-ray-bright stripes ([Eriksen et al v 
2011|) in the South-East region of Tycho (see also |Bykov 
eTaD 20111 In this case, stripes would develop where 
the shock normal is parallel to the large-scale magnetic 
field, and the stripe spacing would correspond to the gy- 
roradius of protons with energy ~ 10 6 GeV, a value com- 
patible with the maximu m energy achieved in Tyc ho as 
inferred from 7-ray data ( |Morlino fc Caprioli||2012| ). 

Though the hybrid simulations presented here cannot 
properly reproduce the large physical scales relevant for 
real SNRs, they qualitatively confirm that the instabil- 
ities driven by the accelerated particles can lead to the 
formation of filaments and knots sharing some proper- 
ties — such as the geometry and the enhancement of the 
magnetic field — consistent with several observational 
signatures detected in young SNRs. 
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